C.............................. PRINTSUBS.FOR ......................
C.. These routines connected to printing FLIP data were moved from 
C.. MONPRN.FOR in July 1998.
C:::::::::::::::::::::::::::::: VPRINT ::::::::::::::::::::::::::::::::
      SUBROUTINE VPRINT(J250,JPRN,IPRIN,IVLPS,MINJ,MAXJ,VBLAT,ISOL,
     >  MAXIT,JBV,ZLO,ZHI,IPP,ALTP,GRD)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER VD,J250
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION IPRIN(55),VPROD(6),VLOSS(6),VCON(VD)
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,D(19),T(2)
     > ,SATX,CHIX,GLONX,GLATX,GRD
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/ND1/NR(2,VD),ZR(VD),TZ(VD),GZ(VD),BZ(VD),VDT,TE(VD),DZ,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      DATA JJJ/0/, VCON/VD*1/

      ID=3
      IF(MINJ.GT.1) ID=4

      !.. Transfer odd N and N2* densities from STR to their correct variables
      CALL TRANF(MINJ,MAXJ,JBV,VCON)

      !..  calculate O(1D), N(2D), NO, and N(4S) if not available in STR .....
      CALL VATMOS(J250,MINJ,MAXJ,GRD,VBLAT,-1,ISOL,MAXIT,IVLPS,JBV)

      IF(I556.NE.-9) WRITE(3,556)
      IF(I556.NE.-9) WRITE(4,556)
      I556=-9
 556  FORMAT(/2X,'**The quantities in this file are given on a vertical'
     > ,1X,'grid at the location specified in the FLIP run file. This'
     > ,/4X,'file is produced when O+, H+, N(2D), NO, N(4S) and'
     > ,1X,'N2(v=0-5) are from time dependent diffusive solutions'/)
C
      DO IS=2,36
        JP=1
        DO J=1,JBV
          IF(ZR(J).GE.ZLO.AND.ZR(J).LE.ZHI) THEN
            IF(JP.EQ.1.OR.MOD(J,IPP).EQ.0) THEN
              IF(IS.EQ.2) CALL PRLN(ID,0,J,JP,3,VPROD,VLOSS,VCON,ICHEM)
              IF(IPRIN(IS).EQ.1) CALL PRLN(ID,IS,J,JP,3,VPROD,VLOSS,VCON
     >        ,ICHEM)
              JP=0
            ENDIF
          ENDIF
        ENDDO
      ENDDO

 202  FORMAT(/6X,'FOR HEATING EFFICIENCY: Thermal electron rates')

      !.. calculate heating efficiency
      IF(IPRIN(35).EQ.1) THEN
        IF(MINJ.EQ.1) THEN
          CALL HEFFIC(ID,JBV,IPP,ZLO,ZHI,MINJ,MAXJ) !.. Northern Hemisphere
        ELSE
          CALL HEFFIC(ID,JBV,IPP,ZLO,ZHI,MINJ,MAXJ) !.. Southern Hemisphere
        ENDIF
      ENDIF
C...... set up parameters for printing N2*
      IF(IPRIN(38).NE.1.AND.IPRIN(39).NE.1) GO TO 485
      IF(IPRIN(38).NE.1) GO TO 463
      CALL PRNN2V(ID,MINJ,MAXJ,-1,4,6,JBV,1,VCON)
      DO 728 J=1,JBV
        IF(ZR(J).LT.ZLO.OR.ZR(J).GT.ZHI) GO TO 728
        IF(MOD(J+IPP-1,IPP).NE.0.AND.J.NE.1) GO TO 728
        CALL PRNN2V(ID,MINJ,MAXJ,J,4,6,JBV,1,VCON)
 728  CONTINUE
C
 463  IF(IPRIN(39).NE.1) GO TO 485
      DO 730 J=1,JBV
      IF(ZR(J).LT.ZLO.OR.ZR(J).GT.ZHI) GO TO 730
      IF(MOD(J+IPP-1,IPP).NE.0.AND.J.NE.1) GO TO 730
      CALL PRVFIJ(ID,MINJ,MAXJ,J,4,6,JBV,1,VCON)
 730  CONTINUE

C
 485  CONTINUE
C........ Vibrationally excited N2
      IF(IPRIN(37).NE.1.OR.IVLPS.LE.1) GO TO 740
      WRITE(ID,733)
 733  FORMAT(/5X,'** temperatures and vibrational populations of N2 **'
     >,3X,'++++ and minor neutral densities +++++'
     >,/3X,'ALT   Tv  Tn   Te   Ti    v=0      v=1      v=2      v=3'
     > ,6X,'v=4      v=5    (N2D)    (N4S)     (NO)     (O1D)')
       DO 735 J=1,JBV
      IF(ZR(J).LT.ZLO.OR.ZR(J).GT.ZHI) GO TO 735
      IF(MOD(J+IPP-1,IPP).NE.0.AND.J.NE.1) GO TO 735
      JCON=J
      IF(MINJ.GT.1) JCON=MAXJ+1-J
      DO 749 I=1,6
 749  D(I)=SNGL(STR(I,JCON))
      CALL TVIB(J,5.0D+2,TV,D)
      CALL AMBS(J,ZR(J),VBLAT,D,T,CHIX,SATX,GLONX,GLATX)
      WRITE(ID,102) ZR(J),NINT(TV),NINT(T(2)),
     >  NINT(TE(J)),NINT(TIZ(J)),
     > (STR(I,JCON),I=1,6),N2D(J),N4S(J),NNO(J),O1D(J)
 735  CONTINUE
 740  CONTINUE


      IF(IPRIN(40).LT.1) RETURN
      DO 377 IS=40,44
        JP=1
      DO 377 J=1,JBV
        IF(ZR(J).LT.ZLO.OR.ZR(J).GT.ZHI) GO TO 377
        IF(MOD(J+IPP-1,IPP).NE.0.AND.J.NE.1) GO TO 377
        !-- store N2* populations for O quenching
        JCON=J
        IF(MINJ.GT.1) JCON=MAXJ+1-J
        DO 44 IL=1,6
           NVS(IL,J)=STR(IL,JCON)
 44     CONTINUE
        IF(IS.EQ.40) CALL PRLN(ID,0,J,JP,3,VPROD,VLOSS,VCON,ICHEM)
        CALL PRLN(ID,IS,J,JP,3,VPROD,VLOSS,VCON,ICHEM)
        JP=0
 377  CONTINUE

 102  FORMAT(F6.1,I5,I5,I5,I5,1P,22E9.2)
      RETURN
      END
C:::::::::::::::::::::::::::::: TVPRINT ::::::::::::::::::::::::::::::::
C... Prints major and minor neutral densities at one altitude
      SUBROUTINE TVPRINT(MINJ,MAXJ,JBV,ALTP,GRD,VBLAT,ISOL,MAXIT,IVLPS)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION VCON(VD)
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,D(19),T(2),SATX,CHIX
      REAL GLONX,GLATX,GRD
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/ND1/NR(2,VD),ZR(VD),TZ(VD),GZ(VD),BZ(VD),VDT,TE(VD),DZ,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)

      DATA JJJ/0/, VCON/VD*1/

      ID=10                 !.. output unit for northern hemisphere
      IF(MINJ.GT.1) ID=11   !.. output unit for southern hemisphere

      !.. Transfer odd N and N2* densities from STR to their correct variables
      CALL TRANF(MINJ,MAXJ,JBV,VCON)

      !..  calculate O(1D), N(2D), NO, and N(4S) if not available in STR .....
      CALL VATMOS(J250,MINJ,MAXJ,GRD,VBLAT,-1,ISOL,MAXIT,IVLPS,JBV)

      !.. Find altitude for printing as a function of time
      DO J=1,JBV-1
        IF(ALTP.LT.(ZR(J+1)+ZR(J))/2.AND.ALTP.GE.(ZR(J)+ZR(J-1))/2) JP=J
      ENDDO

      !.. Take care of ALTP being outside vertical grid range
      IF(ALTP.LE.ZR(1)) JP=1
      IF(ALTP.GE.ZR(JBV)) JP=JBV

      JJJ=JJJ+1
      IF(IVLPS.LE.0) THEN
        IF(JJJ.EQ.1) WRITE(ID,470)
 470    FORMAT(/2X,' *** There are no odd N densities on the input'
     >   ,1X,'FLIP data file ***'/)
        RETURN
      ENDIF

      IF(JJJ.EQ.1) WRITE(10,565) ZR(JP)
      IF(JJJ.EQ.1) WRITE(11,566) ZR(JP)
 565  FORMAT(/9X,'Northern Hemisphere minor and major neutral densities'
     > ,1X,'at alt=',F7.1,' km',
     > /3X,'UT     LT  CHI    N2D      N4S      NO       O        O2'
     > ,7X,'N2      H       He      Tn')
 566  FORMAT(/9X,'Southern Hemisphere minor and major neutral densities'
     > ,1X,'at alt=',F7.1,' km',
     > /3X,'UT     LT  CHI    N2D      N4S      NO       O        O2'
     > ,7X,'N2      H       He      Tn')

      !.. get local time and CHI
      CALL AMBS(JP,ZR(JP),VBLAT,D,T,CHIX,SATX,GLONX,GLATX)

      WRITE(ID,561) SEC/3600.0,SATX,NINT(57.3*CHIX),N2D(JP),
     >  N4S(JP),NNO(JP),ON(JP),O2N(JP),N2N(JP),HN(JP),HE(JP),TZ(JP)
 561  FORMAT(F8.3,F6.2,I4,1P,8E9.2,0P,F7.1,1P,3E9.2)

      RETURN
      END
C::::::::::::::::::::::::::::: BLOCK DATA BLK1  :::::::::::::::::::
      BLOCK DATA BLK1
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      REAL RATFAC
      COMMON/RFAC/RATFAC(99)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
C     change for f90 in Nov, 1999
C      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
C      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      DATA RATFAC/99*1/
      END

C::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE TRANF(MINJ,MAXJ,JB,VCON)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION VCON(VD)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/ND1/NR(2,VD),ZR(VD),TZ(VD),GZ(VD),BZ(VD),VDT,TE(VD),DZ,JTER
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      DO 21 J=1,JB+1

         JCON=J
         IF(MINJ.GT.1) JCON=MAXJ+1-J
         VCON(J)= FITRC(MINJ,MAXJ,Z,RCON,ZR(J))
         DO 25 I=1,4
            NV(I,J)=STR(I+7,JCON)
 25      CONTINUE
         N2D(J)=NV(1,J)
         N4S(J)=NV(2,J)
         NNO(J)=NV(3,J)
         O1D(J)=NV(4,J) * 0.000
 21   CONTINUE
C
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SIMP(IFST,ILST,ISP,H,EN,SUM)
C..... program for integrating an array by simpson's rule written by
C..... P. Richards 12 feb 1980. note that either both ifst @ ilst
C..... are even or they are both odd. h is the distance between grid pts.
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION EN(20,VD)
      EVEN=0.0D0
      ODD=0.0D0
      IF=IFST+1
      IL=ILST-1
C........ add up the odd and even components. fortran indices begin at 1
C........ not 0
      DO 10 J=IF,IL,2
      EVEN=EVEN+EN(ISP,J)
      ODD=ODD+EN(ISP,J+1)
 10   CONTINUE
C......... now form the entire integral 'sum'
      SUM=(EN(ISP,IFST)+4*EVEN+2*ODD+EN(ISP,ILST))*H/3.0
      RETURN
      END
C:::::::::::::::::::::::::::: HEFFIC ::::::::::::::::::::::::::::::::
      SUBROUTINE HEFFIC(ID,JBV,IPP,ZLO,ZHI,MINJ,MAXJ)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER VD,EN_FLAG
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION DELTA(22)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/ENLOSS/ENL(20,VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      DATA SUMHC/0.0D0/,EN_FLAG/0/

      IF(IPP.NE.1) THEN
        WRITE(6,191)
 191    FORMAT(4X,'IPP must be 1 to evaluate heating efficiency'
     > ,/,' integrals')
        RETURN
      ENDIF

      WRITE(ID,203)
 203  FORMAT(/2X,'FOR HEATING EFFICIENCY: Totals of heating and cooling'
     > ,1X,'rates.'
     > /,2X,'O2_DEL = loss due to O2 dissociation;  F_RAD = Forbidden'
     > ,1X,'radiation;   P_RAD = Permitted radiation;  MET_KH = heating'
     > ,1X,'from metastable species'
     > /,2X,'GR_KH = heating from ground state species; O_O_M = 3 body'
     > ,1X,'heating; e_Neut = neutral heating by thermal electrons - '
     > ,1X,'includes electron-ion cooling'  
     > /,2X,'O2_DHT = translational heating from O2 dissociation;' 
     > ,1X,'ION_DEP = total energy deposited in ions;'  
     > ,1X,'S_R = Schumann_Runge energy deposition'
     > /,2X,'N2_DH = heating from dissociation of N2;  e_DEP = total'
     > ,1X,'energy deposited in non-thermal electrons  DIS_DEP ='
     > ,1X,'translational energy of N2 EUV dissociation'
     > /,2X,'EFFIC = (MET_KH + GR_KH + e_Neut + O2_DHT + N2_DH)'
     > ,1X,'/(ION_DEP + S_R + e_DEP + DIS_DEP)'
     > ,/2X,'If the aurora is on, the energy is included. The auroral'
     > ,1X,'and photoelectron energy deposition are in e-DEP.'
     > ,/5X,'ALT',2X,'O2_DEL   F_RAD   P_RAD    MET_KH',4X
     > ,'GR_KH   O_O_M    e_Neut   O2_DHT   ION_DEP    S_R',
     >  4X,'N2_DH     e_DEP  DIS_DEP  Tot_DEP Tot_HEAT  EFFIC')

      DO 97 J=1,JBV
        ENL(14,J)=0.0
        IF(ZR(J).LT.ZLO.OR.ZR(J).GT.ZHI) GO TO 97
        IF(ENL(15,J).GT.1) EN_FLAG=1   !... see if there is any energy deposition

        !.............. heating efficiency .......
        TOTEDEP=ENL(10,J)+ENL(11,J)+ENL(15,J)+ENL(19,J)
        TOTEHEAT=ENL(4,J)+ENL(5,J)+ENL(8,J)+ENL(9,J)+ENL(12,J)

        IF(TOTEHEAT.GT.0.0) ENL(14,J)=TOTEHEAT/TOTEDEP
        IF(ENL(14,J).GT.100.OR.ENL(14,J).LT.0.0001) ENL(14,J)=-99.0

        WRITE(ID,205) ZR(J),(ENL(IK,J),IK=1,5),(ENL(IK,J),IK=7,12)
     >   ,ENL(15,J),ENL(19,J),TOTEDEP,TOTEHEAT,ENL(14,J)
 97   CONTINUE

      DO 85 IS=1,22
 85   DELTA(IS)=0.0
      SUMHC=0.0

      IF(EN_FLAG.NE.1) THEN
        WRITE(6,*) 
     >  ' No Energy input - cannot calculate heating efficiency'
        RETURN
      ENDIF

      !----- integrate assuming exponential parameter variation - convert to ergs
      DO 98 IS=1,20
        DO J=1,JBV-1
          IF(ZR(J).GE.ZLO.AND.ZR(J).LE.ZHI) THEN
            ENPROD=ENL(IS,J)*ENL(IS,J+1)
            IF(ENPROD.GT.0) DELTA(IS)=DELTA(IS)+
     >        DSQRT(ENL(IS,J)*ENL(IS,J+1)) * (ZR(J+1)-ZR(J))*1.6E-7
          ENDIF
        ENDDO
        IF(IS.EQ.4.OR.IS.EQ.5.OR.IS.EQ.8.OR.IS.EQ.9.OR.IS.LT.12)
     >    SUMHC=SUMHC+DELTA(IS)
 98   CONTINUE

      CONSE=0.0D0
      IF(DELTA(10)+DELTA(15)+DELTA(11).GT.1)
     >   CONSE=SUMHC/(DELTA(10)+DELTA(19)+DELTA(15)+DELTA(11))


      !,, WRITE HEATING AND COOLING RATE INTEGRALS

      WRITE(ID,204)
 204  FORMAT(/6X,'FOR HEATING EFFICIENCY: Integrals of the above'
     > ,1X,'heating and cooling rates.'
     > ,/6X,'If an aurora is on, only ~70% of incident energy appears'
     > ,1X,'here because ~30% of auroral energy is backscattered)',//
     > ,5X,'ALT',2X,'O2_DEL   F_RAD   P_RAD    MET_KH',4X
     > ,'GR_KH   O_O_M    e_Neut   O2_DHT   ION_DEP    S_R',
     >  4X,'N2_DH     e_DEP  DIS_DEP  E_Cons')

      WRITE(ID,206) (DELTA(IK),IK=1,5),(DELTA(IK),IK=7,12),DELTA(15),
     >  DELTA(19),CONSE


      DO 333 J=1,JBV
      DO 333 IS=1,20
 333  ENL(IS,J)=0.0D0

 205   FORMAT(F7.1,1P,15E9.2,0P,F8.3)
 206   FORMAT('    ~',1P,14E9.2)
      RETURN
      END
C::::::::::::::::::::: PHEADR::::::::::::::::::::::::::::::::::::::::::::::
C....... Subroutine for printing basic information in each of the
C....... output files
      SUBROUTINE PHEADR(IFILE,PCO,GLONN,GLONF,GLATN,GLATF,ALT,Z0
     >      ,JMAX,SCAL,KVERT)
      DOUBLE PRECISION PCO,SCAL,Z0,GL,ALT
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/SOL/UVFAC(59),EUV
      CHARACTER*60 HEMIS
      IF(IFILE.EQ.3) THEN
         HEMIS = 'NORTHERN hemisphere VERTICAL grid ****'
      ELSE IF(IFILE.EQ.4) THEN
         HEMIS = 'SOUTHERN hemisphere VERTICAL grid ****'
      ELSE IF(IFILE.EQ.9) THEN
         HEMIS = 'NORTHERN hemisphere MAGNETIC field line grid ****'
      ELSE IF(IFILE.EQ.8) THEN
         HEMIS = 'SOUTHERN hemisphere MAGNETIC field line grid ****'
      ELSE
         HEMIS = '****'
      ENDIF
      WRITE(IFILE,112) HEMIS
 112  FORMAT(//5X,'**** Initial Parameters for run with FLIP'
     >   1X,'version 2023-10-13: in the ',A60)
 114  FORMAT(/1X,' YEAR=',I4,',  DAY=',I3,',  L_SHELL=',F8.3,
     >  ',  F10.7=' ,F6.1,',  F10.7A=',F6.1,',  AP=',F6.1
     >  ,' BLON, Z0, JMAX, SCAL = ',3I4,F9.4)
      WRITE(IFILE,114) IDAY/1000,MOD(IDAY,1000),PCO,F107,F107A,AP(1)
     > ,NINT(BLON),NINT(Z0),JMAX,SCAL
 115  FORMAT(/'  Geographic (lat,long):- North (',I3,'N,',I5,'E)'
     > ,2X,',  South (',I3,'S,',I5,'E)  @',I6,'km altitude')
      WRITE(IFILE,115) NINT(GLATN),NINT(GLONN),-NINT(GLATF),NINT(GLONF)
     >  ,NINT(ALT)

      IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3) 
     > WRITE(IFILE,116) UVFAC(9),UVFAC(6),UVFAC(18),UVFAC(33),UVFAC(35)
 116  FORMAT(/'   Using HEUVAC EUV. Initial scaling factors: 303.78='
     >  ,F5.1,', 284.15=',F5.1,  ', 584.33=',F5.1, ', 977.02=',F5.1,
     >  ', 1025.72=',F5.1)
      IF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4)
     > WRITE(IFILE,117) UVFAC(9),UVFAC(6),UVFAC(18),UVFAC(33),UVFAC(35)
 117  FORMAT(/'  Using TIMED-SEE EUV. Initial scaling factors: 303.78='
     >  ,F5.1,', 284.15=',F5.1,  ', 584.33=',F5.1, ', 977.02=',F5.1,
     >  ', 1025.72=',F5.1)
      IF(NINT(UVFAC(58)).GE.0) 
     > WRITE(IFILE,118) UVFAC(9),UVFAC(6),UVFAC(18),UVFAC(33),UVFAC(35)
 118  FORMAT(/'   User provided EUV scaling factors: 303.78='
     >  ,F5.1,', 284.15=',F5.1,  ', 584.33=',F5.1, ', 977.02=',F5.1,
     >  ', 1025.72=',F5.1)

      !... IF(KVERT.LE.0) WRITE(IFILE,120) 
 120  FORMAT(///'  <<< There are no altitude profiles on the vertical'
     > ,1X,'grid because the continuity and momentum equations'
     > ,/6X,'for odd N were not solved on the FLIP run.'
     > ,//1X,' >>> Look in the field line files for odd N chemical'
     > ,1X,'equilibrium solutions instead.'///) 

      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
        SUBROUTINE PRNN2V(ID,MINJ,MAXJ,J,ILJ,IPR,JB,IC,VCON)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION B(22),VCON(VD),VPROD(6),VLOSS(6),C(22)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      IF(J.LE.0) WRITE(ID,95)
 95   FORMAT(/11X,'Total production and loss of N2 vibrational quanta'
     > ,/22X,'Production',27X,':',9X,'Loss'
     > ,/3X,'ALT    e-therm    O(1D)    NO+N4S     N2(A)   e > 1eV'
     > ,4X,'N2+(v)    N2(v)+O   N2(v)+O+')
      IF(J.LE.0) RETURN
C
      DO 10 IK=1,22
10    C(IK)=0.0D0
C
      JCON=J
      IF(MINJ.GT.1) JCON=MAXJ+1-J
C
      DO 13 IS=1,4
         NV(IS,J)=STR(IS+7,JCON)
 13   CONTINUE
C
      N2D(J)=STR(8,JCON)
      N4S(J)=STR(9,JCON)
      NNO(J)=STR(10,JCON)
      O1D(J)=STR(11,JCON) * 0.0000
      IF(O1D(J).LE.0) CALL PRLN(1,0,J,0,3,VPROD,VLOSS,VCON,0)
C
      DO 15 IS=1,IPR
         NV(IS,J)=STR(IS,JCON)
 15   CONTINUE
C
      DO 20 IS=2,IPR
      CALL PRLO(MINJ,MAXJ,ILJ,IS,J,IPR,PROD,LOSS,B)
C
      C(1)=C(1)+(IS-1)*(B(1)-B(2))
      C(2)=C(2)+(IS-1)*(B(8)-B(7))
      IF(IS.EQ.3) C(3)=C(3)+(IS-1)*B(9)
      IF(IS.EQ.5) C(4)=C(4)+(IS-1)*B(9)
      IF(IS.EQ.IPR) C(5)=C(5)+(IS-1)*B(13)
      C(6)=C(6)+(IS-1)*B(12)
      C(7)=C(7)+(IS-1)*B(14)
C....      C(8)=C(8)+(IS-1)*B(13?)
 20   CONTINUE
C
      WRITE(ID,99) ZR(J),C(1),C(3),C(4),C(5),C(7),C(8),C(2),C(6)
 99   FORMAT(F7.1,1P,19E10.2)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE HTPRIN(JJJ,JPRIN,IFILE,UTHRS,SATN,CHI,HFLUX1,HFLUX2
     > ,HIFLX1,HIFLX2)
      DOUBLE PRECISION HFLUX1,HFLUX2,HIFLX1,HIFLX2
C     change UTHRS from double to real for f90, Nov, 1999
      REAL UTHRS
      DOUBLE PRECISION TEFLUX,TFLUX,AVEE,TEEHT,FRPAS
      COMMON/EFLUX/TEFLUX,TFLUX,AVEE,TEEHT,FRPAS
      IF(JJJ.EQ.1.OR.JPRIN.GT.0) WRITE(IFILE,365)
 365  FORMAT(/9X,'These parameters are concerned with the flow of heat
     > to and from the plasmasphere(POS) from both ends of the'
     > ,/5X,'flux tube. HTFLUX=total thermal electron heat flux from
     > POS. AVEE=average energy of escaping photoelectrons'
     > ,/5X,'TEFLUX=total energy flux carried into the POS by
     > photoelectrons. TFLUX=the total particle flux carried'
     > ,/5X,'into the POS. TEEHT=total amount of electron-electron
     > heating in the POS.'
     > ,//4X,'UT',4X,'LT',3X,'CHI   AVEE   HTFLUX   TEFLUX   TFLUX'
     > ,' TEEHT')
      WRITE(IFILE,361) UTHRS,SATN,NINT(57.3*CHI),AVEE,HFLUX1+HFLUX2
     > ,TEFLUX,TFLUX,TEEHT,HIFLX1+HIFLX2
 361  FORMAT(F7.2,F6.2,I4,F7.1,1P,22E10.2)
      RETURN
      END
C:::::::::::::::::::::::::: FNDTIM :::::::::::::::::::::::::::::::::::::::
C..... This subroutine reads the F8.DAT data file and matches the
C..... given time (TSTART) with the nearest time on the file then stores
C..... the old value in TSAVE before returning the nearest time in TSTART
C..... and TSTOP. If TSTART is not in the range between the first file
C..... time and the last file time the program stops
C..... Program was completely rewritten for simplicity in December 2008 by
C..... P. Richards. This was achieved by storing the UT in an array
      SUBROUTINE FNDTIM(TSTART,TSTOP,TSAVE,IDAY,PCO,BLON)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER IDIM,NHRS
      PARAMETER (IDIM=999999)
      REAL SEC,F107,F107A,AP(7),BLON,EUV,GLONN,GLONF,PLAT,PLON, 
     > GLATN,GLATF,GRD,UVFAC(59),RATFAC(99),TIMEON,TIMEOF,AURFAC
     > ,TSTART,TSTOP,TSAVE
      REAL UTHRS(IDIM)
      DOUBLE PRECISION Q(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)

      IFLAG=0   !.. Set FLAG to indicate end of file

      !.. Set Trap invalid start and stop time
      IF(TSTART.GT.TSTOP) THEN
        WRITE(6,'(/A)') 
     >   ' ** ERROR: Requested TSTART > TSTOP in FLIPRINT'
        STOP
      ENDIF

      !.. save stop and start times
      STARTSAVE=TSTART
      STOPSAVE=TSTOP

      CALL READ_Q(SCALK, Q)  !.. get field line coordinate system

      !..... Main loop to read UT into an array
      DO J=1,IDIM
        CALL RF8DAT(IVLPS,JQ,IW,PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN
     >     ,GLATF,ZLBDY,FPAS,UV,HFAC,JMAX,IDAY,I1,I2,F107,F107A,AP
     >     ,BLON,EUV,Z,UVFAC,RATFAC,TIMEON,TIMEOF,AURFAC,PLAT,PLON)
        CALL VERGRD(IVLPS,GRD) !.. vertical grid to initialize JBV,JVB,JNB
        CALL RF8D2(JTX,JK,JC,NMF2N,NMF2F,NOPEQ,NHPEQ,TIEQ,TEEQ
     >     ,UNN,UNC,IVLPS,IFLAG,JQ,SEC,N,JMAX,TI,U,STR,RCON,AP
     >     ,NHEPEQ,U,F107,NNPEQ)
        IF(IFLAG.EQ.1) GO TO 20
        UTHRS(J)=SEC/3600.
      ENDDO

 20   NHRS=J-1   !.. # of UTs
      REWIND 1   !.. Rewind unit 1 for MONPRN 
      IFLAG=0
       
      !..  Test if TSTART is too small
      IF(TSTART.LT.UTHRS(2)) THEN
        TSTART=UTHRS(2)
        STARTSAVE=TSTART
        WRITE(6,'(/A,F10.2,A)') 
     >    ' ** TSTART < 1st UT time on file, reset to',TSTART,' UT'
      ENDIF


      !..  Test if TSTART is too Big
      IF(STARTSAVE.GT.UTHRS(NHRS)) THEN
        !.. See if last time is close enough to requested start time
        IF(STARTSAVE-UTHRS(NHRS).LE.0.25) THEN
          STARTSAVE=UTHRS(NHRS)
        ELSE 
          WRITE(6,*) ' ** ERROR: TSTART > last time on file'
          STOP
        ENDIF
      ENDIF

      !..  Test if TSTOP is too small
      IF(TSTOP.LT.TSTART) THEN
        WRITE(6,'(/A)') ' ** ERROR: TSTOP < TSTART'
        STOP
      ENDIF

      !..  Test if TSTOP is too Big
      IF(TSTOP.GT.UTHRS(NHRS)) THEN
        !.. See if last time is close enough to requested start time
        IF(STOPSAVE-UTHRS(NHRS).LE.0.25) THEN
          STOPSAVE=UTHRS(NHRS)
        ELSE 
          WRITE(6,'(/A,F10.2,A)') 
     >     ' ** ERROR: Requested TSTOP > last UT time on file',
     >      UTHRS(NHRS),' UT'
          STOP
        ENDIF
      ENDIF
      TSTART=UTHRS(NHRS)
      TSTOP=UTHRS(NHRS)
      !.. Look for the closest UT for TSTART and TSTOP
      DO J=1,NHRS-1
        IF(STARTSAVE.LT.(UTHRS(J+1)+UTHRS(J))/2.AND.
     >     STARTSAVE.GE.(UTHRS(J)+UTHRS(J-1))/2) TSTART=UTHRS(J)
        IF(STOPSAVE.LT.(UTHRS(J+1)+UTHRS(J))/2.AND. 
     >     STOPSAVE.GE.(UTHRS(J)+UTHRS(J-1))/2)  TSTOP=UTHRS(J)
      ENDDO

      RETURN
      END
C::::::::::::::::::::::::::: FITRC ::::::::::::::::::::::::::::::::::
      FUNCTION FITRC(MINJ,MAXJ,ZP,T,Z)
C--- This program was written by P. Richards in August 1991
C...... Linear interpolation of vibrational N2 factor (y = A.x + B)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION ZP(401),T(401)
C
C--- Use bisection to find the nearest altitude ZP(K) to desired
C--- altitude (Z)
      CALL BISPLT(MINJ,MAXJ,Z,ZP,K)
C
C....... now do linear fit
 2    IF(K.EQ.MAXJ) THEN
         FITRC=T(K)
      ELSE IF(T(K).LT.0.OR.T(K+1).LT.0) THEN
         FITRC=0.0
      ELSE
         FACTOR=(Z-ZP(K))/(ZP(K+1)-ZP(K))
         FITRC=T(K)+FACTOR*(T(K+1)-T(K))
      ENDIF
 4     RETURN
       END
C:::::::::::::::::::::::::::: TRSTR :::::::::::::::::::::::::::::
C---- This subroutine transfers the vertical grid odd N densities to
C---- the field line grid for use in CMINOR. Written by P. Richards
C---- April 1994
      SUBROUTINE TRSTR(IVLPS)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER IVERT,VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION ZRD(401)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/VODDN/VN2D(401),VN4S(401),VNNO(401),VO1D(401),IVERT
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/ND1/NR(2,VD),ZR(VD),TZ(VD),GZ(VD),BZ(VD),VDT,TE(VD),DZ,JTER
      COMMON/VAGRID/JBV,JVB,JNB
C
      IVERT=IVLPS
      IF(IVERT.LE.0) RETURN
C--- A dummy vertical altitude grid is created so that values in STR can
C--- be interpolated
      DO 10 J=1,JBV       
         ZRD(J)=ZR(J)
         ZRD(JMAX+1-J)=ZR(J)
 10   CONTINUE
C
C---- Now extract odd N densities from STR and interpolate
      DO 21 J=1,JMAX/2
         IF(Z(J).LT.ZR(JBV)) THEN
            VN2D(J)=FITSTR(8,1,JBV,ZRD,STR,Z(J))
            VN4S(J)=FITSTR(9,1,JBV,ZRD,STR,Z(J))
            VNNO(J)=FITSTR(10,1,JBV,ZRD,STR,Z(J))
C..            VO1D(J)=FITSTR(11,1,JBV,ZRD,STR,Z(J))
            JCON=JMAX+1-J
            VN2D(JCON)=FITSTR(8,JMAX+1-JBV,JMAX,ZRD,STR,Z(JCON))
            VN4S(JCON)=FITSTR(9,JMAX+1-JBV,JMAX,ZRD,STR,Z(JCON))
            VNNO(JCON)=FITSTR(10,JMAX+1-JBV,JMAX,ZRD,STR,Z(JCON))
C..            VO1D(JCON)=FITSTR(11,JMAX+1-JBV,JMAX,ZRD,STR,Z(JCON))
C...           WRITE(68,90) J, NINT(Z(J)), VN2D(J),VN4S(J),VNNO(J)
           JC=JCON
C...           WRITE(69,90) JC, NINT(Z(JC)), VN2D(JC),VN4S(JC),VNNO(JC)
 90        FORMAT(1X,I4,I7,1P,9E10.2)  
       ENDIF
 21   CONTINUE
C
      RETURN
      END
C::::::::::::::::::::::::::: FITSTR ::::::::::::::::::::::::::::::::::
      FUNCTION FITSTR(I,JMIN0,JMAX0,ZP,T,Z)
C...... Linear interpolation for temperatures (y = A.x + B)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION ZP(401), T(12,401)
C
C--- Use bisection to find the nearest altitude ZP(K) to desired 
C--- altitude (Z)
      CALL BISPLT(JMIN0,JMAX0,Z,ZP,I1)
C
C....... now do linear fit
 2    IF(I1.EQ.JMAX0) THEN
         FITSTR=T(I,I1)
      ELSE IF(T(I,I1).LT.0.OR.T(I,I1+1).LT.0) THEN
         FITSTR=0.0
      ELSE
         FACTOR=(Z-ZP(I1))/(ZP(I1+1)-ZP(I1))
         FITSTR=T(I,I1)+FACTOR*(T(I,I1+1)-T(I,I1))
      ENDIF
 4     RETURN
       END

C:::::::::::::::::::: GET_AE_OPLUS :::::::::::::::::::::::::::::
C.... This routine is used to get O+ data from a file to
C.... replace the FLIP values in the chemical scheme
C.... It takes the FLIP altitude and returns an interpolated
C.... value. Bad or no data are returned as 0.
      SUBROUTINE GET_AE_DENS(FLIP_ALT,
     >                          UTHRS,   !... UT in hours
     >                          GLATD,   !... latitude (deg)
     >                          AE_OP,   !... AE O+ density
     >                          AE_OX,   !... AE O density
     >                          AE_N2,   !... AE N2 density
     >                          AE_O2,   !... AE O2 density
     >                          AE_NO,   !... AE N2 density
     >                          AE_TE,   !... AE electron temp
     >                          AE_TI)   !... AE ion temp
      IMPLICIT NONE
      INTEGER ALTDIM   !.... Array dimensions
      INTEGER MAXALT   !.... Index of maximum altitude in array
      INTEGER JCALL    !.... counts number of calls
      INTEGER DIMALT
      PARAMETER (ALTDIM=99)
      REAL UTHRS,TSTOP,TSTART !... universal time in hours
      REAL GLATD,LATDAT       !... Flip and data latitude
      REAL AE_ALT(ALTDIM)     !... Altitudes of AE data
      REAL AE_OPLUS(ALTDIM)   !... AE O+ density
      REAL AE_ONEUT(ALTDIM)   !... AE O density
      REAL AE_N2NEUT(ALTDIM)  !... AE N2 density
      REAL AE_O2NEUT(ALTDIM)  !... AE N2 density
      REAL AE_NONEUT(ALTDIM)  !... AE NO density
      REAL AE_TEMPTE(ALTDIM)  !... AE electron temperature
      REAL AE_TEMPTI(ALTDIM)  !... AE ion temperature
      REAL FTAE               !... Interpolating function
      REAL FLIP_ALT           !... FLIP altitude for interpolation
      REAL OPLUS              !... Temporary O+ variable, error check
      !... AE densities and temperatures to return
      REAL AE_OP, AE_OX, AE_N2, AE_O2, AE_NO, AE_TE, AE_TI 
      DATA JCALL/0/
      DIMALT=ALTDIM

      !... if this is the first time, read the AE data file
      JCALL=JCALL+1
      IF(JCALL.EQ.1) 
     >   CALL READ_AE_DATA(DIMALT,AE_ALT,AE_OPLUS,AE_ONEUT
     >       ,AE_N2NEUT,AE_O2NEUT,AE_NONEUT,AE_TEMPTE,AE_TEMPTI,MAXALT,
     >       TSTART,TSTOP,LATDAT)

      AE_OP=0.0
      AE_OX=0.0
      AE_N2=0.0
      AE_O2=0.0
      AE_NO=0.0
      AE_TE=0.0
      AE_TI=0.0

      IF(LATDAT*GLATD.LT.0) RETURN    !.. check correct hemisphere
      IF(UTHRS.LT.TSTART.OR.UTHRS.GT.TSTOP) RETURN
     
      !... get interpolated value of O+ if in AE alt range
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_OP= FTAE(1,ALTDIM,AE_ALT,AE_OPLUS,FLIP_ALT)
      !... get interpolated value of O+ if in AE alt range
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_OX= FTAE(1,ALTDIM,AE_ALT,AE_ONEUT,FLIP_ALT)
      !... get interpolated value of O+ if in AE alt range
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_N2= FTAE(1,ALTDIM,AE_ALT,AE_N2NEUT,FLIP_ALT)
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_O2= FTAE(1,ALTDIM,AE_ALT,AE_O2NEUT,FLIP_ALT)
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_NO= FTAE(1,ALTDIM,AE_ALT,AE_NONEUT,FLIP_ALT)
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_TE= FTAE(1,ALTDIM,AE_ALT,AE_TEMPTE,FLIP_ALT)
      IF(FLIP_ALT.GT.AE_ALT(1).AND.FLIP_ALT.LE.AE_ALT(MAXALT))
     >   AE_TI= FTAE(1,ALTDIM,AE_ALT,AE_TEMPTI,FLIP_ALT)

      !... Only return the AE value if it is good, else 0.
      IF(AE_OP.LT.100.OR.AE_OP.GT.5.0E6)  AE_OP=0.0
      IF(AE_OX.LT.100.OR.AE_OX.GT.1.0E20) AE_OX=0.0
      IF(AE_N2.LT.100.OR.AE_N2.GT.1.0E20) AE_N2=0.0
      IF(AE_O2.LT.100.OR.AE_O2.GT.1.0E20) AE_O2=0.0
      IF(AE_NO.LT.10.OR.AE_NO.GT.1.0E20)  AE_NO=0.0
      IF(AE_TE.LT.300.OR.AE_TE.GT.1.0E4) AE_TE=0.0
      IF(AE_TI.LT.300.OR.AE_TI.GT.1.0E4) AE_TI=0.0
      RETURN
      END

C:::::::::::::::::::: READ_AE_DATA :::::::::::::::::::::::::::
C... This subroutine actually reads the AE data from the file
C... Data must either increase monotonically ordecrease monotonically
      SUBROUTINE READ_AE_DATA(ALTDIM,   !... Array dimension in
     >                       AE_ALT,    !... AE altitude array out
     >                       AE_OPLUS,  !... AE O+ density
     >                       AE_ONEUT,  !... AE O density
     >                       AE_N2NEUT, !... AE N2 density
     >                       AE_O2NEUT, !... AE O2 density
     >                       AE_NONEUT, !... AE NO density
     >                       AE_TEMPTE, !... AE NO density
     >                       AE_TEMPTI, !... AE NO density
     >                       MAXALT,    !... Index of max alt
     >                       TSTART,    !... time to start using data
     >                       TSTOP,     !... time to stop using data
     >                       LATDAT)    !... latitude of data
      IMPLICIT NONE
      INTEGER INCOLS         !... dimension for # columns to read
      INTEGER ALTDIM         !... Array dimensions in GET_AE_DENS
      INTEGER DIMALT         !... Array dimensions in READ_AE
      PARAMETER (DIMALT=99,INCOLS=11)
      INTEGER I,J,JREV       !... loop control variable
      INTEGER MAXALT         !... Index of maximum altitude in array
      INTEGER FEXISTS        !... Indicates if file exists
      INTEGER NCOLS          !... Number of columns to read
      REAL TSTART,TSTOP      !... start and stop times for data
      REAL LATDAT            !... latitude in data
      REAL COLS(INCOLS)      !... Column switches to read and multiplier
      REAL TALT(DIMALT),AE_ALT(DIMALT)    !... Altitudes of AE data
      REAL TOPD(DIMALT),AE_OPLUS(DIMALT)  !... AE O+ density
      REAL TOXD(DIMALT),AE_ONEUT(DIMALT)  !... AE O density
      REAL TN2D(DIMALT),AE_N2NEUT(DIMALT) !... AE N2 density
      REAL TO2D(DIMALT),AE_O2NEUT(DIMALT) !... AE N2 density
      REAL TNOD(DIMALT),AE_NONEUT(DIMALT) !... AE NO density
      REAL TTE(DIMALT),AE_TEMPTE(DIMALT)  !... AE NO density
      REAL TTI(DIMALT),AE_TEMPTI(DIMALT)  !... AE NO density
      REAL AEDATA(INCOLS)    !... Temp array to read a line of data. 
      INTEGER K  !&&&&&&&&&&
      MAXALT=0
      !... Using DIMALT in here because ALTDIM generates a large
      !... executable file for some reason
      IF(DIMALT.NE.ALTDIM) THEN
         WRITE(6,*) ' ** ERROR: wrong dimension for DIMALT in READ_AE'
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF

      !... see if file exists and signal with negative TSTOP
      CALL TFILE(47,FEXISTS)
      TSTOP=-1000
      IF(FEXISTS.EQ.0) RETURN

      !... remove headers from file, reading to column info
 10   READ(47,*,ERR=10,END=91) TSTART,TSTOP

      !... Error trap for empty file - Xiang Zhao on May 09,2002
      GO TO 92      !... if success skip over this block
 91   TSTOP=-1000  !.. return TSTOP as flag
      RETURN

 92   WRITE(6,877)
 877  FORMAT(/3X,' Data are being read from file to replace FLIP'
     >  ,1X,'O+, O, and N2 densities'/)

      READ(47,*) LATDAT
      READ(47,*) (COLS(J),J=1,INCOLS)
      WRITE(6,*) ' TSTART    TSTOP    LATDAT  (COLS(J),J=1,INCOLS)'
      WRITE(6,'(3F9.2,11F6.2)') TSTART,TSTOP,LATDAT,(COLS(J),J=1,INCOLS)
 
      !... Determine number of columns of data to read
      DO J=1,INCOLS
         IF(COLS(J).GT.0.05) NCOLS=J
      ENDDO
       WRITE(6,*) ' **** AEAEAE ****'
      WRITE(6,*) NCOLS
      !... Read the AE data into temporary arrays
      DO J=1,999
         READ(47,*,ERR=95,END=93) (AEDATA(I), I=1,NCOLS)
         TOXD(J)=0.0
         TN2D(J)=0.0
         TO2D(J)=0.0
         TOPD(J)=0.0
         TNOD(J)=0.0
         TTE(J)=0.0
         TTI(J)=0.0
         IF(COLS(1).GT.0.05) TALT(J)=AEDATA(1)
         IF(COLS(2).GT.0.05) TOXD(J)=COLS(2)*AEDATA(2)
         IF(COLS(3).GT.0.05) TN2D(J)=COLS(3)*AEDATA(3)
         IF(COLS(4).GT.0.05) TO2D(J)=COLS(4)*AEDATA(4)
         IF(COLS(5).GT.0.05) TOPD(J)=COLS(5)*AEDATA(5)
         IF(COLS(6).GT.0.05) TNOD(J)=COLS(6)*AEDATA(6)
         IF(COLS(7).GT.0.05) TTE(J)=COLS(7)*AEDATA(7)
         IF(COLS(8).GT.0.05) TTI(J)=COLS(8)*AEDATA(8)
      ENDDO

 93   MAXALT=J-1

      !.. If altitudes not already in increasing, reverse order, then
      !.. test for monotonic increasing
      DO J=1,MAXALT
         JREV=J
         IF(TALT(1).GE.TALT(MAXALT)) JREV=MAXALT+1-J
         AE_ALT(J)=TALT(JREV)
         AE_OPLUS(J)=TOPD(JREV)
         AE_ONEUT(J)=TOXD(JREV)
         AE_N2NEUT(J)=TN2D(JREV)
         AE_O2NEUT(J)=TO2D(JREV)
         AE_NONEUT(J)=TNOD(JREV)
         AE_TEMPTE(J)=TTE(JREV)
         AE_TEMPTI(J)=TTI(JREV)
         !.. check for monotonically increasing order
         IF(J.GT.1.AND.AE_ALT(J).LE.AE_ALT(J-1)) THEN
            WRITE(6,191) 
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF
      ENDDO

      CLOSE(47)
      RETURN
 95   WRITE(6,*) ' Error in data file'
      CALL RUN_ERROR    !.. print ERROR warning in output files
      STOP

 191  FORMAT(/' ** ERROR in Altitude data file unit = 47'
     > ,1X,'Altitudes must be monotonically increasing or decreasing')
      END
C:::::::::::::::::::::: FUNCTION FTAE ::::::::::::::::::::::::::::::::
C...... Simple linear interpolating program for AE data
      REAL FUNCTION FTAE(BEGI,   !... smallest array index
     >                   IDIM,   !... largest array index
     >                     ZP,   !... Altitude array
     >                      T,   !... array of valuse to interp
     >                      Z)   !... alt for interpolated value
      IMPLICIT NONE
      INTEGER BEGI,IDIM,I
      REAL ZP(IDIM),T(IDIM),Z,FACTOR
      DO 1 I=BEGI,IDIM-1 
       IF(Z.GE.ZP(I).AND.Z.LT.ZP(I+1)) GOTO 2
 1     CONTINUE
 2     FACTOR=(Z-ZP(I))/(ZP(I+1)-ZP(I))
       FTAE=T(I)+FACTOR*(T(I+1)-T(I))
       IF (Z .EQ. ZP(IDIM)) FTAE = T(IDIM)
 4    RETURN
      END
